function R2 = IPCA_R2(Gamma,F,R,Z,IsF,X,W,varargin)

% Gamma is L*(K+M)
% F is K*T
% R is N*T
% Z is N*T*L
% IsF is N*T
% X is L*T
% W is L*L*T
% PSF is M*T

% Standard IPCA call as IPCA_R2(Gamma,F,R,Z,IsF,X,W)
% IPCA with PSF call as IPCA_R2(Gamma,F,R,Z,IsF,X,W,PSF) - Gamma should be L*(K+M)
% IPCA with alpha call as IPCA_R2(Gamma,F,R,Z,IsF,X,W,ones(1,T)) - Gamma should be L*(K+1)
% IPCA with ONLY PSF call as IPCA_R2(Gamma,[],R,Z,IsF,X,W,PSF) - Gamma should be L*M

if nargin>7 && ~isempty(varargin{1})
    PSF = varargin{1};
    PSF_version = true;
    M = size(PSF,1);
    [L, Ktilde] = size(Gamma);
    K = Ktilde - M;
    if K>0
        GammaBeta = Gamma(:,1:K);
        LambdaF = mean(F,2,'omitnan');
        SigmaF = cov(F','omitrows');
        w = SigmaF\LambdaF; w = w/sum(w);
        Tangency_F = w'*F;
        TangencySR_F = sqrt(12)*mean(Tangency_F,'omitnan')/std(Tangency_F,'omitnan');
    end
    GammaDelta = Gamma(:,K+1:end);
    LambdaPSF = mean(PSF,2,'omitnan');
    SigmaPSF = cov(PSF','omitrows');
    w = SigmaPSF\LambdaPSF; w = w/sum(w);
    if mean(w'*PSF,2,'omitnan') < 0
        w = -w;
    end
    Tangency_PSF = w'*PSF;
    TangencySR_PSF = sqrt(12)*mean(Tangency_PSF,'omitnan')/std(Tangency_PSF,'omitnan');
else
    PSF_version = false;
    [L, K] = size(Gamma);
    GammaBeta = Gamma(:,1:K);
    LambdaF = mean(F,2,'omitnan');
    SigmaF = cov(F','omitrows');
    w = SigmaF\LambdaF; w = w/sum(w);
    Tangency_F = w'*F;
    TangencySR_F = sqrt(12)*mean(Tangency_F,'omitnan')/std(Tangency_F,'omitnan');
end
[N, T] = size(R);

%% R2 for individiual assets
numer1 = NaN(N,T); numer2 = NaN(N,T); denom = NaN(N,T);
for t = 1:T-1
    Zt = squeeze(Z(:,t,:));
    Rt1 = R(:,t+1);
    ValidIndex = IsF(:,t+1);

    Exp1 = zeros(N,1); Exp2 = zeros(N,1);
    if K>0
        Betat = Zt*GammaBeta;
        Ft1 = F(:,t+1);
        Exp1 = Exp1 + Betat*Ft1;
        Exp2 = Exp2 + Betat*LambdaF;
    end
    if PSF_version
        Deltat = Zt*GammaDelta;
        PSFt1 = PSF(:,t+1);
        Exp1 = Exp1 + Deltat*PSFt1;
        Exp2 = Exp2 + Deltat*LambdaPSF;
    end
    
    Err1t1 = Rt1 - Exp1;
    Err2t1 = Rt1 - Exp2;
    
    numer1(ValidIndex,t+1) = Err1t1(ValidIndex).^2;
    numer2(ValidIndex,t+1) = Err2t1(ValidIndex).^2;
    denom(ValidIndex,t+1) = Rt1(ValidIndex).^2;
end
TotalR2_i = 1 - sum(numer1(:),'omitnan')/sum(denom(:),'omitnan');
PredR2_i = 1 - sum(numer2(:),'omitnan')/sum(denom(:),'omitnan');
R2i = 1 - sum(numer1,2,'omitnan')./sum(denom,2,'omitnan'); R2i(~isfinite(R2i)|sum(~isnan(numer1),2)<6) = NaN; Ti = sum(isfinite(numer1),2); TimeSeriesR2_i = sum(Ti.*R2i,'omitnan')/sum(Ti,'omitnan');
R2x = 1 - sum(numer1,1,'omitnan')./sum(denom,1,'omitnan'); R2x(~isfinite(R2x)) = NaN; CrossSectionR2_i = mean(R2x,'omitnan');

%% Relative pricing error for individiual assets
numer1 = NaN(N,T); denom = NaN(N,T);
for t = 1:T-1
    Zt = squeeze(Z(:,t,:));
    Rt1 = R(:,t+1);
    ValidIndex = IsF(:,t+1);

    Exp1 = zeros(N,1);
    if K>0
        Betat = Zt*GammaBeta;
        Ft1 = F(:,t+1);
        Exp1 = Exp1 + Betat*Ft1;
    end
    if PSF_version
        Deltat = Zt*GammaDelta;
        PSFt1 = PSF(:,t+1);
        Exp1 = Exp1 + Deltat*PSFt1;
    end
    Err1t1 = Rt1 - Exp1;
    
    numer1(ValidIndex,t+1) = Err1t1(ValidIndex);
    denom(ValidIndex,t+1) = Rt1(ValidIndex);
end
Alpha = mean(numer1,2,'omitnan');
Rbar = mean(denom,2,'omitnan');
RelPrcErr_i = sum(Alpha.^2,'omitnan')/sum(Rbar.^2,'omitnan');

%% R2 for managed portfolios
numer1 = NaN(L,T); numer2 = NaN(L,T); denom = NaN(L,T);
for t = 1:T-1
    Wt = W(:,:,t);
    Xt1 = X(:,t+1);

    Exp1 = zeros(L,1); Exp2 = zeros(L,1);
    if K>0
        Ft1 = F(:,t+1);
        Exp1 = Exp1 + Wt*GammaBeta*Ft1;
        Exp2 = Exp2 + Wt*GammaBeta*LambdaF;
    end
    if PSF_version
        PSFt1 = PSF(:,t+1);
        Exp1 = Exp1 + Wt*GammaDelta*PSFt1;
        Exp2 = Exp2 + Wt*GammaDelta*LambdaPSF;
    end
    
    Err1t1 = Xt1 - Exp1;
    Err2t1 = Xt1 - Exp2;
    
    numer1(:,t+1) = Err1t1.^2;
    numer2(:,t+1) = Err2t1.^2;
    denom(:,t+1) = Xt1.^2;
end
%numer1(end,:) = []; numer2(end,:) = []; denom(end,:) = []; % the last 'managed' portfolio is just the average return; remove it
TotalR2_x = 1 - sum(numer1(:),'omitnan')/sum(denom(:),'omitnan');
PredR2_x = 1 - sum(numer2(:),'omitnan')/sum(denom(:),'omitnan');
R2i = 1 - sum(numer1,2,'omitnan')./sum(denom,2,'omitnan'); Ti = sum(isfinite(numer1),2); TimeSeriesR2_x = sum(Ti.*R2i,'omitnan')/sum(Ti,'omitnan');
R2x = 1 - sum(numer1,1,'omitnan')./sum(denom,1,'omitnan'); CrossSectionR2_x = mean(R2x,'omitnan');

%% Relative pricing error for managed portfolios
numer1 = NaN(L,T); denom = NaN(L,T);
for t = 1:T-1
    Wt = W(:,:,t);
    Xt1 = X(:,t+1);

    Exp1 = zeros(L,1);
    if K>0
        Ft1 = F(:,t+1);
        Exp1 = Exp1 + Wt*GammaBeta*Ft1;
    end
    if PSF_version
        PSFt1 = PSF(:,t+1);
        Exp1 = Exp1 + Wt*GammaDelta*PSFt1;
    end
    Err1t1 = Xt1 - Exp1;
    
    numer1(:,t+1) = Err1t1;
    denom(:,t+1) = Xt1;
end
%numer1(end,:) = []; denom(end,:) = []; % the last 'managed' portfolio is just the average return; remove it
Alpha = mean(numer1,2,'omitnan');
Rbar = mean(denom,2,'omitnan');
RelPrcErr_x = sum(Alpha.^2,'omitnan')/sum(Rbar.^2,'omitnan');


%--------------------------------------------------------------------------
R2.TotalR2_i = TotalR2_i;
R2.PredR2_i = PredR2_i;
R2.TotalR2_x = TotalR2_x;
R2.PredR2_x = PredR2_x;
R2.TimeSeriesR2_i = TimeSeriesR2_i; 
R2.CrossSectionR2_i = CrossSectionR2_i;
R2.TimeSeriesR2_x = TimeSeriesR2_x;
R2.CrossSectionR2_x = CrossSectionR2_x;
R2.RelPrcErr_i = RelPrcErr_i;
R2.RelPrcErr_x = RelPrcErr_x;
if K>0
    R2.TangencySR_F = TangencySR_F;
end
if PSF_version
    R2.TangencySR_PSF = TangencySR_PSF;
end

return